{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "0b9ffbb2-8d33-eb0c-5a42-0733bdc3b74b",
    "_uuid": "e14f65723f5ac826104e861c45e58525740d8b2a"
   },
   "source": [
    "# Bike Sharing 数据集上的回归分析\n",
    "\n",
    "1、 任务描述 请在Capital Bikeshare （美国Washington, D.C.的一个共享单车公司）提供的自行车数据上进行回归分析。训练数据为2011年的数据，要求预测2012年每天的单车共享数量。\n",
    "\n",
    "原始数据集地址：http://archive.ics.uci.edu/ml/datasets/Bike+Sharing+Dataset 1) 文件说明 \n",
    "day.csv: 按天计的单车共享次数（作业只需使用该文件） \n",
    "hour.csv: 按小时计的单车共享次数（无需理会） \n",
    "readme：数据说明文件\n",
    "\n",
    "2) 字段说明 \n",
    "Instant记录号 \n",
    "Dteday：日期 \n",
    "Season：季节（1=春天、2=夏天、3=秋天、4=冬天） \n",
    "yr：年份，(0: 2011, 1:2012) \n",
    "mnth：月份( 1 to 12) \n",
    "hr：小时 (0 to 23) （只在hour.csv有，作业忽略此字段） \n",
    "holiday：是否是节假日 \n",
    "weekday：星期中的哪天，取值为0～6 \n",
    "workingday：是否工作日 1=工作日 （是否为工作日，1为工作日，0为非周末或节假日 weathersit：天气（1：晴天，多云 ",
    "2：雾天，阴天 ",
    "3：小雪，小雨 ",
    "4：大雨，大雪，大雾） temp：气温摄氏度 \n",
    "atemp：体感温度 \n",
    "hum：湿度 \n",
    "windspeed：风速 \n",
    "\n",
    "casual：非注册用户个数 \n",
    "registered：注册用户个数 \n",
    "cnt：给定日期（天）时间（每小时）总租车人数，响应变量y\n",
    "casual、registered和cnt三个特征均为要预测的y，作业里只需对cnt进行预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "_cell_guid": "2ba3154a-c2aa-3158-1984-63ad2c0c786a",
    "_uuid": "5eb696b95780825e94ddb49787f9fa339fc3833b"
   },
   "outputs": [],
   "source": [
    "# 导入必要的工具包\n",
    "# 数据读取及基本处理\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "#导入模型\n",
    "from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV\n",
    "\n",
    "#评价回归预测模型的性能\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "#可视化\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "_cell_guid": "21fa35be-878b-b4f2-ef6e-68dc070b8bfa",
    "_uuid": "73aee228226be55c0c8a6e4fcbf818c56cd94926",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season_1</th>\n",
       "      <th>season_2</th>\n",
       "      <th>season_3</th>\n",
       "      <th>season_4</th>\n",
       "      <th>mnth_1</th>\n",
       "      <th>mnth_2</th>\n",
       "      <th>mnth_3</th>\n",
       "      <th>mnth_4</th>\n",
       "      <th>mnth_5</th>\n",
       "      <th>...</th>\n",
       "      <th>weekday_5</th>\n",
       "      <th>weekday_6</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>holiday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>yr</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.355170</td>\n",
       "      <td>0.373517</td>\n",
       "      <td>0.828620</td>\n",
       "      <td>0.284606</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.379232</td>\n",
       "      <td>0.360541</td>\n",
       "      <td>0.715771</td>\n",
       "      <td>0.466215</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.171000</td>\n",
       "      <td>0.144830</td>\n",
       "      <td>0.449638</td>\n",
       "      <td>0.465740</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.175530</td>\n",
       "      <td>0.174649</td>\n",
       "      <td>0.607131</td>\n",
       "      <td>0.284297</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.209120</td>\n",
       "      <td>0.197158</td>\n",
       "      <td>0.449313</td>\n",
       "      <td>0.339143</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 35 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant  season_1  season_2  season_3  season_4  mnth_1  mnth_2  mnth_3  \\\n",
       "0        1         1         0         0         0       1       0       0   \n",
       "1        2         1         0         0         0       1       0       0   \n",
       "2        3         1         0         0         0       1       0       0   \n",
       "3        4         1         0         0         0       1       0       0   \n",
       "4        5         1         0         0         0       1       0       0   \n",
       "\n",
       "   mnth_4  mnth_5  ...   weekday_5  weekday_6      temp     atemp       hum  \\\n",
       "0       0       0  ...           0          1  0.355170  0.373517  0.828620   \n",
       "1       0       0  ...           0          0  0.379232  0.360541  0.715771   \n",
       "2       0       0  ...           0          0  0.171000  0.144830  0.449638   \n",
       "3       0       0  ...           0          0  0.175530  0.174649  0.607131   \n",
       "4       0       0  ...           0          0  0.209120  0.197158  0.449313   \n",
       "\n",
       "   windspeed  holiday  workingday  yr   cnt  \n",
       "0   0.284606        0           0   0   985  \n",
       "1   0.466215        0           0   0   801  \n",
       "2   0.465740        0           1   0  1349  \n",
       "3   0.284297        0           1   0  1562  \n",
       "4   0.339143        0           1   0  1600  \n",
       "\n",
       "[5 rows x 35 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 读入数据\n",
    "data = pd.read_csv(\"FE_day.csv\")\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "a1ae46d1-8787-67a3-2f69-6497c97320eb",
    "_uuid": "16fa6600f9877c9607060f9e50de6eaa9bc1c766"
   },
   "source": [
    "数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#得到输出y\n",
    "y = data['cnt']   \n",
    "X = data.drop(['cnt'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train samples: (584, 34)\n"
     ]
    }
   ],
   "source": [
    "# 分割训练数据和测试数据\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2,random_state = 0)\n",
    "\n",
    "print(\"train samples:\" ,X_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jhony/anaconda3/lib/python3.7/site-packages/pandas/core/frame.py:3697: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  errors=errors)\n"
     ]
    }
   ],
   "source": [
    "#ID列不参与预测\n",
    "X_train.drop(['instant'], axis=1, inplace = True)\n",
    "X_test.drop(['instant'], axis=1, inplace = True)\n",
    "\n",
    "#保存特征名字\n",
    "feat_names = X_train.columns\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "9fe1c63e-3803-feac-362c-631399fdb8ec",
    "_uuid": "431749cca6d8fbbc2d1d3732baeb6d63590c6e2e"
   },
   "source": [
    "1. 最小二乘线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "_cell_guid": "101cf15a-9006-ac9a-8abe-756371fd8b1a",
    "_uuid": "38562579d6306877a067189f25ce15ececbe99d2",
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE on Training set : 756.1481528838043\n",
      "RMSE on Test set : 785.9189127479793\n",
      "r2_score on Training set : 0.8420556108191108\n",
      "r2_score on Test set : 0.8547736196975189\n"
     ]
    }
   ],
   "source": [
    "# 1. 生成学习器实例\n",
    "lr = LinearRegression()\n",
    "\n",
    "#2. 在训练集上训练学习器\n",
    "lr.fit(X_train, y_train)\n",
    "\n",
    "#3.用训练好的学习器对训练集/测试集进行预测\n",
    "y_train_pred = lr.predict(X_train)\n",
    "y_test_pred = lr.predict(X_test)\n",
    "\n",
    "\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\", rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\", r2_score_train)\n",
    "print(\"r2_score on Test set :\", r2_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "season_1       -8.183624e+15\n",
       "season_2       -8.183624e+15\n",
       "season_3       -8.183624e+15\n",
       "season_4       -8.183624e+15\n",
       "mnth_1          1.883302e+16\n",
       "mnth_2          1.883302e+16\n",
       "mnth_3          1.883302e+16\n",
       "mnth_4          1.883302e+16\n",
       "mnth_5          1.883302e+16\n",
       "mnth_6          1.883302e+16\n",
       "mnth_7          1.883302e+16\n",
       "mnth_8          1.883302e+16\n",
       "mnth_9          1.883302e+16\n",
       "mnth_10         1.883302e+16\n",
       "mnth_11         1.883302e+16\n",
       "mnth_12         1.883302e+16\n",
       "weathersit_1    2.618741e+16\n",
       "weathersit_2    2.618741e+16\n",
       "weathersit_3    2.618741e+16\n",
       "weekday_0      -2.929307e+16\n",
       "weekday_1       2.095560e+16\n",
       "weekday_2       2.095560e+16\n",
       "weekday_3       2.095560e+16\n",
       "weekday_4       2.095560e+16\n",
       "weekday_5       2.095560e+16\n",
       "weekday_6      -2.929307e+16\n",
       "temp            2.882000e+03\n",
       "atemp           1.108500e+03\n",
       "hum            -1.842000e+03\n",
       "windspeed      -1.450250e+03\n",
       "holiday        -5.024867e+16\n",
       "workingday     -5.024867e+16\n",
       "yr              1.933312e+03\n",
       "dtype: float64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "系数的值非常大。由于特征之间强相关，OLS模型的性能并不好"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "a70803d8-638c-bdfa-6897-aed8ce78f475",
    "_uuid": "2938aee78a9f61d6ca59345d6fa0793e2422f34c"
   },
   "source": [
    "2.岭回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 1\n",
      "RMSE on Training set : 754.0366623762052\n",
      "RMSE on Test set : 776.9753607134023\n",
      "r2_score on Training set : 0.8429364764067757\n",
      "r2_score on Test set : 0.8580600896877546\n"
     ]
    }
   ],
   "source": [
    "# 1. 设置超参数搜索范围，生成学习器实例\n",
    "# RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, normalize=False, scoring=None, cv=None, gcv_mode=None, store_cv_values=False)\n",
    "alphas = [0.01, 0.1, 1, 10, 100, 1000]\n",
    "ridge = RidgeCV(alphas = alphas,cv = 5)\n",
    "\n",
    "# 2. 用训练数据对模型进行训练\n",
    "ridge.fit(X_train, y_train)\n",
    "\n",
    "#通过交叉验证得到的最佳超参数alpha\n",
    "alpha = ridge.alpha_\n",
    "print(\"Best alpha :\", alpha)\n",
    "\n",
    "#模型评估\n",
    "y_train_pred = ridge.predict(X_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "\n",
    "y_test_pred = ridge.predict(X_test)\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "\n",
    "print(\"RMSE on Training set :\", rmse_train)\n",
    "print(\"RMSE on Test set :\" ,rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\" ,r2_score_train)\n",
    "print(\"r2_score on Test set :\" ,r2_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "season_1        -816.946054\n",
       "season_2          61.439057\n",
       "season_3          97.124573\n",
       "season_4         658.382424\n",
       "mnth_1          -384.604861\n",
       "mnth_2          -255.398769\n",
       "mnth_3           185.553801\n",
       "mnth_4            -6.012570\n",
       "mnth_5           437.523845\n",
       "mnth_6           104.170972\n",
       "mnth_7          -307.981162\n",
       "mnth_8           129.351248\n",
       "mnth_9           670.800620\n",
       "mnth_10          235.474439\n",
       "mnth_11         -428.094512\n",
       "mnth_12         -380.783051\n",
       "weathersit_1     757.284463\n",
       "weathersit_2     341.945472\n",
       "weathersit_3   -1099.229935\n",
       "weekday_0       -166.919676\n",
       "weekday_1       -155.850725\n",
       "weekday_2        -48.116950\n",
       "weekday_3         -4.759450\n",
       "weekday_4         64.295894\n",
       "weekday_5         87.446625\n",
       "weekday_6        223.904283\n",
       "temp            1892.310062\n",
       "atemp           1605.349484\n",
       "hum            -1520.738791\n",
       "windspeed      -1336.567102\n",
       "holiday         -202.590794\n",
       "workingday       145.606188\n",
       "yr              1952.328786\n",
       "dtype: float64"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "相比OLS，岭回归模型增加了L2正则，系数值进行了收缩。\n",
    "由于增加正则限制了模型复杂度，相比比OLS模型，岭回归模型在训练集上和测试集上的误差有所减小。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "2da33461-fb1a-c09a-b0fe-d20e8c8f52e3",
    "_uuid": "d2e8bab153344d6ecf2e48909b76bd8f387e2710"
   },
   "source": [
    "3. Lasso模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "_cell_guid": "8525724a-fb77-66d4-e06f-f3365fbd8ec3",
    "_uuid": "f8f41944e405693ba867158353ffedcd9ecd4a07",
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best alpha : 2.0321278108385803\n",
      "cv of rmse : 809.677633728043\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEKCAYAAAC7c+rvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt8VdWd9/HPL1fuhEtQDCCgiILKpRGxVsdbFa0tbdVW26nUsbXTqc+0nXaqbWdqLzPzajtTfWovdmz1qdZWy6O9UCuPpYpap4qGOxGBAAIhgYQQciGQ6+/546zgISbkEDnZ55x836/X9uyz9tp7/TYHz4+19jp7m7sjIiIShayoAxARkYFLSUhERCKjJCQiIpFREhIRkcgoCYmISGSUhEREJDJKQiIiEhklIRERiYySkIiIRCYn6gBS3dixY33y5MlRhyEiklZWrly5z90Le6unJNSLyZMnU1JSEnUYIiJpxcx2JFJPw3EiIhIZJSEREYmMkpCIiERGSUhERCKjJCQiIpFREhIRkcgoCYmISGSUhERE5C2+/+ctvLS1JuntKAmJiMhR9tQd5n8/s5lX39if9LaUhERE5ChLN1TiDtecMz7pbSkJiYjIUZ5aX8mZJw/n9HHDkt6WkpCIiByxp+4wJTtq+6UXBElMQmY2yMxeMbO1ZlZqZt8I5T83s+1mtiYss0O5mdm9ZlZmZuvMbG7csRaZ2ZawLIorf4eZrQ/73GtmFspHm9myUH+ZmY3qrQ0REenfoThIbk+oGbjM3WcBs4EFZjY/bPtnd58dljWh7GpgWlhuA+6DWEIB7gLOB+YBd3UmlVDntrj9FoTyO4Fn3H0a8Ex432MbIiIS89T6Sqaf1D9DcZDEJOQxjeFtblj8GLssBB4O+70MFJjZeOAqYJm773f3WmAZsYQ2Hhjh7i+5uwMPA++PO9ZDYf2hLuXdtSEiMuDtre/foThI8jUhM8s2szVAFbFEsiJs+vcwHHaPmeWHsiJgV9zu5aHsWOXl3ZQDnOTulQDhdVwvbYiIDHhL18eG4t5z7sn91mZSk5C7t7v7bGACMM/Mzga+DJwJnAeMBu4I1a27Q/Sh/FgS2sfMbjOzEjMrqa6u7uWQIiKZ4Y9HhuKG91ub/TI7zt0PAM8BC9y9MgyHNQP/h9h1Hoj1SibG7TYBqOilfEI35QB7O4fZwmtVL210jfd+dy929+LCwl6fTisikvb6e1Zcp2TOjis0s4KwPhi4Ang9LjkYsWs1G8IuS4Cbwwy2+UBdGEp7GrjSzEaFCQlXAk+HbQ1mNj8c62bg93HH6pxFt6hLeXdtiIgMaE+uq8Ad3jurf5NQThKPPR54yMyyiSW7xe7+pJk9a2aFxIbG1gB/H+o/BVwDlAFNwC0A7r7fzL4FvBrqfdPdO+8l8Wng58BgYGlYAL4NLDazW4GdwA3HakNEZKBbsraCs4tGMLWwf2bFdUpaEnL3dcCcbsov66G+A5/pYduDwIPdlJcAZ3dTXgNcfjxtiIgMVNv3HWRdeR1fveasfm9bd0wQERnglqypwAyu7eehOFASEhEZ0NydJWt3c97k0YwfObjf21cSEhEZwF6rrGdr9UEWzj4lkvaVhEREBrAlayvIyTKuOTuam8coCYmIDFAdHc4f1lRw0bSxjBqaF0kMSkIiIgPU6l0HqKg7zHtnRTMUB0pCIiID1tL1leRmG1fMOCmyGJSEREQGIHdn6YY9XDStkBGDciOLQ0lIRGQAWr+7jt0HDrHg7P67Y3Z3lIRERAagpRv2kJNlXBnhUBwoCYmIDDjuztL1lVxw2hgKhkQzK66TkpCIyADz+p4G3qhp4uqIfhsUT0lIRGSAWbq+kiyDK2dGOxQHSkIiIgPO0g17mDdlNGOH5UcdipKQiMhAsn3fQbZUNabEUBwoCYmIDChrdx0A4PypoyOOJEZJSERkACmtqCMvJ4vT+vkJqj1REhIRGUBKK+o58+Th5Ganxtd/akQhIiJJ5+6UVtQz85QRUYdyhJKQiMgAsfvAIeoOtTLjlJFRh3KEkpCIyABRWlEPoJ6QiIj0v9KKerIMzjpZSUhERPrZaxV1TC0cxuC87KhDOUJJSERkgEi1SQmgJCQiMiDsP9hCZd3hgZOEzGyQmb1iZmvNrNTMvhHKp5jZCjPbYma/NrO8UJ4f3peF7ZPjjvXlUL7JzK6KK18QysrM7M648uNuQ0Qkk5VW1AEwM4VmxkFye0LNwGXuPguYDSwws/nAd4B73H0aUAvcGurfCtS6++nAPaEeZjYDuBGYCSwAfmxm2WaWDfwIuBqYAdwU6nK8bYiIZLpUnBkHSUxCHtMY3uaGxYHLgMdD+UPA+8P6wvCesP1yM7NQ/pi7N7v7dqAMmBeWMnff5u4twGPAwrDP8bYhIpLRSivqKSoYHPlD7LpK6jWh0GNZA1QBy4CtwAF3bwtVyoGisF4E7AII2+uAMfHlXfbpqXxMH9roGvdtZlZiZiXV1dV9O3kRkRTyWkUdM1KsFwRJTkLu3u7us4EJxHouZ3VXLbx21yPxE1h+rDaOLnC/392L3b24sLCwm11ERNJHU0sb2/YdTLmhOOin2XHufgB4DpgPFJhZTtg0AagI6+XARICwfSSwP768yz49le/rQxsiIhlrY2UD7qk3KQGSOzuu0MwKwvpg4ApgI7AcuD5UWwT8PqwvCe8J2591dw/lN4aZbVOAacArwKvAtDATLo/Y5IUlYZ/jbUNEJGOtCc8QOndC6iWhnN6r9Nl44KEwiy0LWOzuT5rZa8BjZvZvwGrggVD/AeAXZlZGrHdyI4C7l5rZYuA1oA34jLu3A5jZ7cDTQDbwoLuXhmPdcTxtiIhkstU7aykqGMxJIwZFHcpbmDoCx1ZcXOwlJSVRhyEi0mcXfvtZZk8q4EcfmdtvbZrZSncv7q2e7pggIpLB9tYfZveBQ8ydNCrqULqlJCQiksFW74xdD5ozqSDiSLqnJCQiksFW76wlLzsrJadng5KQiEhGW73zADOLRpCfkzqPb4inJCQikqFa2ztYt/sAcyam5vUgUBISEclYr1c2cLi1I2WvB4GSkIhIxlq9qxZI3UkJoCQkIpKxVu88wLjh+RQVDI46lB4pCYmIZKjVO2uZM6mAVH5ijZKQiEgGqmls5o2aJuak6I9UOykJiYhkoCM/Up2YuteDQElIRCQjrdheQ152FrOUhEREpL+9tK2GOZMKGJSbmj9S7aQkJCKSYeoOtVJaUc8Fp42JOpReKQmJiGSYV7bvxx3mT1USEhGRfvbS1hryc7JS+keqnZSEREQyzMvbanjHqaNS9qal8ZSEREQyyIGmFjbuqU+LoThQEhIRySgvb4tdD0qHSQmgJCQiklFe3lbD4NxsZk1I/etBoCQkIpJRXt5WQ/HkUeTlpMfXe3pEKSIivappbOb1PQ1pcz0IlIRERDLGiu37gfT4fVAnJSERkQzxwuZqhufncO6EkVGHkrCkJSEzm2hmy81so5mVmtlnQ/nXzWy3ma0JyzVx+3zZzMrMbJOZXRVXviCUlZnZnXHlU8xshZltMbNfm1leKM8P78vC9sm9tSEiks7cneWbqrj4jEJys9Onf5HMSNuAL7j7WcB84DNmNiNsu8fdZ4flKYCw7UZgJrAA+LGZZZtZNvAj4GpgBnBT3HG+E441DagFbg3ltwK17n46cE+o12MbyfsjEBHpH69V1rO3vplLphdGHcpxSVoScvdKd18V1huAjUDRMXZZCDzm7s3uvh0oA+aFpczdt7l7C/AYsNBijwq8DHg87P8Q8P64Yz0U1h8HLg/1e2pDRCStLX+9CoBLpo+LOJLj0y99tjAcNgdYEYpuN7N1ZvagmXU+9q8I2BW3W3ko66l8DHDA3du6lB91rLC9LtTv6VgiImlt+aZqzp0wksLh+VGHclySnoTMbBjwBPA5d68H7gNOA2YDlcD3Oqt2s7v3obwvx+oa821mVmJmJdXV1d3sIiKSOmoPtrB6Zy2XplkvCJKchMwsl1gC+qW7/wbA3fe6e7u7dwA/5c3hsHJgYtzuE4CKY5TvAwrMLKdL+VHHCttHAvuPcayjuPv97l7s7sWFhek1vioiA88LW6rpcLj0TCWhI8I1mAeAje5+d1z5+LhqHwA2hPUlwI1hZtsUYBrwCvAqMC3MhMsjNrFgibs7sBy4Puy/CPh93LEWhfXrgWdD/Z7aEBFJW8tfr2LM0DzOLUqfqdmdcnqv0mcXAh8D1pvZmlD2FWKz22YTGwZ7A/gUgLuXmtli4DViM+s+4+7tAGZ2O/A0kA086O6l4Xh3AI+Z2b8Bq4klPcLrL8ysjFgP6Mbe2hARSUftHc7zm6u59MxxZGV1d8UhtVmsgyA9KS4u9pKSkqjDEBHp1sodtVx331/5wU1zeO+sU6IO5wgzW+nuxb3VS59fNImIyFs8s3Ev2VnGxdPS8/q1kpCISJpyd55aX8k7TxvDyCG5UYfTJ0pCIiJpqrSinjdqmnjPOeN7r5yilIRERNLUH9dXkp1lXDnz5KhD6TMlIRGRNBQ/FDd6aF7U4fSZkpCISBoqrahnR5oPxYGSkIhIWuocirsqjYfiQElIRCTtuDt/XBcbihuVxkNxoCQkIpJ2NuyuZ+f+Jq49N72H4kBJSEQk7Ty5riI2K25Geg/FgZKQiEhaaW3v4IlVu7l0+ri0H4qD40hCZvYuM7slrBeGu1CLiEg/em5TNfsam/lQ8YSoQzkhEkpCZnYXsTtWfzkU5QKPJCsoERHp3uKSXYwdlp+Wzw7qTqI9oQ8A7wMOArh7BTA8WUGJiMhbVTUc5tnXq7hubhG52ZlxNSXRs2gJD4VzADMbmryQRESkO79dtZv2DueG4om9V04TiSahxWb238Qep/1J4M/EHs0tIiL9wN1ZXLKL4lNHcfq4YVGHc8Ik9GRVd/8vM3s3UA9MB77m7suSGpmIiByxamctW6sP8t3rTos6lBMqoSQUht+edfdlZjYdmG5mue7emtzwREQE4NFXdjEkL5v3ZMAPVOMlOhz3ApBvZkXEhuJuAX6erKBERORN+xqbWbK2gg/OLWJofkJ9h7SRaBIyd28CPgj8wN0/AMxIXlgiItLpkZd30NLWwS0XZt7PMxNOQmZ2AfBR4I+hLLPSsYhICjrc2s4jL+/gsjPHcVph5kxI6JRoEvoscCfwG3cvDXdLeDZ5YYmICMCStRXsa2zh7zKwFwSJ92aagA7gJjP7W8AIvxkSEZHkcHcefHE7Z548nAtPHxN1OEmRaBL6JfBFYAOxZCQiIkn21601vL6nge9edy5mFnU4SZHocFy1u//B3be7+47O5Vg7mNlEM1tuZhvNrNTMPhvKR5vZMjPbEl5HhXIzs3vNrMzM1pnZ3LhjLQr1t5jZorjyd5jZ+rDPvRY+pb60ISKSan76l22MGZrH+2afEnUoSZNoErrLzH5mZjeZ2Qc7l172aQO+4O5nAfOBz5jZDGLXlp5x92nAM+E9wNXAtLDcBtwHsYQC3AWcD8wLsYwK+9wX6nbutyCUH1cbIiKpZs2uAzy3qZq/e9cUBuVmRx1O0iSahG4BZhP7kn9vWK491g7uXunuq8J6A7ARKAIWAg+Fag8B7w/rC4GHPeZlYrcIGg9cBSxz9/3uXgssAxaEbSPc/aVwX7uHuxzreNoQEUkp3//zZkYNyWXROydHHUpSJXpNaJa7n9PXRsxsMjAHWAGc5O6VEEtUZtZ5P/IiYFfcbuWh7Fjl5d2U04c2Kvt6biIiJ9rqnbUs31TNlxZMZ1iG/Ti1q0R7Qi+HobTjZmbDgCeAz7l7/bGqdlPmfSg/ZjiJ7GNmt5lZiZmVVFdX93JIEZET6/vPbGHUkFxuvmBy1KEkXaJJ6F3AGjPbFC7orzezdb3tZGa5xBLQL939N6F4b+cQWHitCuXlQPz9yScAFb2UT+imvC9tHMXd73f3YncvLiws7O00RUROmNU7a3luUzWfvHhqxveCIPEktIDYxfwrefN60HuPtUOYqfYAsNHd747btATonOG2CPh9XPnNYQbbfKAuDKk9DVxpZqPChIQrgafDtgYzmx/aurnLsY6nDRGRlHD3snAtaAD0giDxRzkcczp2Dy4EPgasN7M1oewrwLeJPZ/oVmAncEPY9hRwDVBG7Mext4S295vZt4BXQ71vuvv+sP5pYjdSHQwsDQvH24aISCp4blMVf9myj395z1kZd6PSnlhsYpn0pLi42EtKSqIOQ0QyXFt7B1d//y+0tHfwp89fTH5Oek/LNrOV7l7cW73MeEi5iEiae+zVXWypauTLV5+V9gnoeCgJiYhErP5wK/cs28z5U0Zz1cyTog6nXykJiYhE7EfLy9jf1MK/XjsjY+8R1xMlIRGRCG3a08CDL27nurkTOLtoZNTh9DslIRGRiLR3OHc8sY7hg3L5yjVnRR1OJJSEREQi8ouX3mDNrgN87doZjB6aF3U4kVASEhGJwO4Dh/ju05v4mzMKWZjBj2rojZKQiEg/c3f+5bfrAfj3D5w94CYjxFMSEhHpZ0vWVrB8UzVfvHI6E0YNiTqcSCkJiYj0o9qDLXzzD68xa2JBxj8rKBED4+ZEIiIp4t+f2kjdoVYe+eA5ZGcN3GG4TuoJiYj0kxe37OPxleXcdvFUzho/IupwUoKSkIhIPzjU0s5Xf7eeyWOG8I+XT4s6nJSh4TgRkX5wz583s6OmiUc/OZ9BuQPnBqW9UU9IRCTJ1u46wM/+so2PnD+JC04bE3U4KUVJSEQkiVraOvjS4+sYN3wQd159ZtThpBwNx4mIJNF9z21l094GHlhUzIhBuVGHk3LUExIRSZKNlfX8cPkWFs4+hcvPGljPCUqUkpCISBK0tHXwT4vXMnJwHne9d2bU4aQsDceJiCTB95/ZzMbKen52c/GAvUN2ItQTEhE5wVbtrOW+57ZywzsmcMUMDcMdi5KQiMgJdKilnS8uXsv4kYP52ntnRB1OytNwnIjICfStP77G9pqD/PLW8xmu2XC9Uk9IROQE+VPpHn61Yie3XTSVd54+Nupw0oKSkIjICbC3/jB3PLGOmaeM4AtXTo86nLSRtCRkZg+aWZWZbYgr+7qZ7TazNWG5Jm7bl82szMw2mdlVceULQlmZmd0ZVz7FzFaY2RYz+7WZ5YXy/PC+LGyf3FsbIiJvR0eH84XFaznU2s73b5xDXo7+fZ+oZP5J/RxY0E35Pe4+OyxPAZjZDOBGYGbY58dmlm1m2cCPgKuBGcBNoS7Ad8KxpgG1wK2h/Fag1t1PB+4J9Xps4wSfs4gMQI+s2MGLZfv42rUzOX3csKjDSStJS0Lu/gKwP8HqC4HH3L3Z3bcDZcC8sJS5+zZ3bwEeAxZa7IHslwGPh/0fAt4fd6yHwvrjwOWhfk9tiIj0WXltE99Z+joXTRvLTfMmRh1O2omiz3i7ma0Lw3WjQlkRsCuuTnko66l8DHDA3du6lB91rLC9LtTv6VhvYWa3mVmJmZVUV1f37SxFJOO5O1/57QYc+I8PnEPs37tyPPo7Cd0HnAbMBiqB74Xy7j4570N5X4711kL3+9292N2LCwsLu6siIsITq3bzwuZq7lhwJhNHD4k6nLTUr0nI3fe6e7u7dwA/5c3hsHIgvh87Aag4Rvk+oMDMcrqUH3WssH0ksWHBno4lInLcqhoO860nX6P41FF8bP6pUYeTtvo1CZnZ+Li3HwA6Z84tAW4MM9umANOAV4BXgWlhJlwesYkFS9zdgeXA9WH/RcDv4461KKxfDzwb6vfUhojIcfv6klIOtbbz7evOJStLw3B9lbQ7JpjZo8AlwFgzKwfuAi4xs9nEhsHeAD4F4O6lZrYYeA1oAz7j7u3hOLcDTwPZwIPuXhqauAN4zMz+DVgNPBDKHwB+YWZlxHpAN/bWhojI8Xi6dA9Prd/DP181XbPh3iaLdRKkJ8XFxV5SUhJ1GCKSIuoOtfLuu59nzLB8ltx+IbnZ+k1Qd8xspbsX91ZPf3oiIsfh20s3sq+xme9cd44S0AmgP0ERkQS9tLWGR1/ZxScumsq5EwqiDicjKAmJiCTgUEs7d/5mHZNGD+HzV5wRdTgZQ49yEBFJwPf+tIkdNU08+sn5DM7THb9OFPWERER6sWpnLQ/8z3Y+ev4kLjhtTNThZBQlIRGRY2hua+dLj69j/IhB3Hn1mVGHk3E0HCcicgzf+9Nmyqoa+fkt5+lJqUmgnpCISA+Wv17F/S9s46PnT+KS6eOiDicjKQmJiHSjsu4Q/7R4DWeePJx/vXZG7ztInygJiYh00dbewWcfW0NzWwc//MhcBuVqNlyy6JqQiEgX//WnzbyyfT93f2iW7g2XZOoJiYjEWfzqLn7y/FY+cv4kPjh3QtThZDwlIRGR4K9l+/jKb9dz0bSxfON9M6MOZ0BQEhIRAcqqGvn7R1YyZexQfviRubo5aT/Rn7KIDHg7a5r42AMryM3O4sGPn8fIwfo9UH/RxAQRGdDKa5u46acvc6i1nV99Yj4TRw+JOqQBRT0hERmwKusOcdNPX6bhcCuP3Ho+M04ZEXVIA456QiIyIJVVNbLowVeoP9TKI584n7OLRkYd0oCkJCQiA07JG/v5xMMl5GRl8atPzuecCUpAUVESEpEBZen6Sj736zUUFQzmob+bp2tAEVMSEpEBoa29g//80yb++/ltzJ1UwM8WncfooXlRhzXgKQmJSMarajjM//rValZs38/fzp/Ev147g/wc3Q8uFSgJiUhGe7p0D1/97Xoam9u4+0OzdCueFKMkJCIZqa6pla//oZTfrt7NzFNG8MsPzWb6ycOjDku6SNrvhMzsQTOrMrMNcWWjzWyZmW0Jr6NCuZnZvWZWZmbrzGxu3D6LQv0tZrYorvwdZrY+7HOvmVlf2xCRzNHR4fxmVTlX3PM8S9ZW8NnLp/G7z1yoBJSikvlj1Z8DC7qU3Qk84+7TgGfCe4CrgWlhuQ24D2IJBbgLOB+YB9zVmVRCndvi9lvQlzZEJHOsL6/j+p/8lX9avJZTCgbzu3+4kM+/+wzdBy6FJW04zt1fMLPJXYoXApeE9YeA54A7QvnD7u7Ay2ZWYGbjQ91l7r4fwMyWAQvM7DlghLu/FMofBt4PLD3eNty98kSet4j0v817G7hn2WaWbtjD2GF5fPf6c7l+7gSysizq0KQX/X1N6KTOL313rzSzzoe2FwG74uqVh7JjlZd3U96XNpSERNLU+vI67v/LNp5cV8HQvBz+8bLT+cTFUxkxSDcgTRepMjGhu3+ueB/K+9LGWyua3UZsyI5Jkyb1clgR6U9t7R08+3oVD7y4nRXb9zM0L5tPXXwan7p4KqP0u5+0099JaG/nEFgYbqsK5eXAxLh6E4CKUH5Jl/LnQvmEbur3pY23cPf7gfsBiouLe0tuItIPtu87yOKSXTyxspyqhmaKCgbz1WvO4sPzJqrnk8b6OwktARYB3w6vv48rv93MHiM2CaEuJJGngf+Im4xwJfBld99vZg1mNh9YAdwM/KAvbSTxXEXkbdpZ08Qf11fy1PpK1u+uI8vg0unj+NB5E7n8zHHkaMJB2ktaEjKzR4n1YsaaWTmxWW7fBhab2a3ATuCGUP0p4BqgDGgCbgEIyeZbwKuh3jc7JykAnyY2A28wsQkJS0P5cbUhIqmjrb2DVTsP8Mzre1n+ehWb9zYCMGtiAV+++kzeP6eIk0YMijhKOZEsNllMelJcXOwlJSVRhyGSsQ40tfCXLft4ZuNenttczYGmVnKyjPOnjubS6eNYcPbJTBilm4ymGzNb6e7FvdVLlYkJIjJAHGhqYfXOA7y8vYa/ltWwoaIOdxg1JJfLzhzHFWedxEXTxjJc13kGBCUhEUma5rZ2SivqWV9ex7ryOtaWH6CsKjbElpttzJk4is9dfgbvmjaG2RNHka3f9Qw4SkIicsIcbm1nXXkdK7bV8NK2GlbuqKW5rQOAscPyOHdCAR+YU8TcSaOYNXEkQ/L0FTTQ6W+AiPRZS1sHK3fU8mJZNSu27WddeR0t7bGkM2P8CP52/qmcN3k0syaO5OQRgwi3eBQ5QklIRBLW2t4R6+lsr2HFtv28sn0/h1rbyc4yzikayccvnMx5k0dTfOoo/XBUEqIkJCLHdKilnec3V7F0wx6e2VhFY3MbANPGDeOG4glcNK2Q+VNHayKB9ImSkIi8RVXDYZa/XsWy16p4sayaw60djBqSy7Xnjudvzihk3pTRjBmWH3WYkgGUhESE6oZmVu+s5a9ba3hpaw2b9jYAUFQwmA8VT2TBzJOZN2W07lAgJ5ySkMgA09HhbKlq5JXtNbz6Ri2rd9Wya/8hAAblZnHe5NEsnHMKl04fx5knD9dkAkkqJSGRDNdwuJX1u+tYvfMAK3fUsmpnLQeaWgE4aUQ+7zh1FDfPn8ycSQWcM2Ek+TnZEUcsA4mSkEiG6OhwymsPsWlvA5v3NrCxsp7Sinq27zt4pM7p44Zx5YyTOG/yaM6fMoaJowerpyORUhISSTMNh1vZUdPEtn0HeWPfQbZWN1JW1cjW6kYOt3YcqTdh1GDOPmUk180tYmbRSOZMLKBgiKZNS2pREhJJQfWHW3lj30G27zvI1uqDbKtuZOf+JnbubzoylNbplJGDmHbScOZPHcO0ccOYfvJwpp00nGH5+t9bUp/+lopE5FBLOzv2H+SNfU1s29fI9uqDR3o3NQdbjtTLMpgwaginjhnCe84Zz8TRQ5g0eghTC4dy6uihDM7TNRxJX0pCIknQ3uHUNDZT1dBMZd1h9tQdorLuMLsPHKK89hDltU3srW8+ap/C4flMGTuUd884iSljhzJ57FCmjB3KqWOGaLKAZCwlIZEEtLV3UH+4jQNNLdQ2tVB7sJX9B1uoOdjC/oPN1BxsYV9jC9UNzexrbKamsZmOLo/qyskyxhcMoqhgMO86vZDJY4Zw6tihTB4zhMljh+oR1TIgKQklSUtbB00tbRgGBhZ7IcssrMdeIbYtywwD7MgrmrXUDXenw2M9jfYOp92dtvYO2sL71vYOWttjZS1hvaWtg5a2Dg63ttMcXg+3tXOopZ3Dre00tXQubRxsbqexuY3G5jYaDrfScLiNhsNtR25V051BuVmMGZrP2OFntiwRAAAHx0lEQVT5FBUMYtaEkYwbnk/hiEEUDstn/MhBjB85iDHD8vWoApEulISS5E+v7eH2X61+28eJT1BHElhcEnszqUFWlh2dzCx2PcGw2Gvc/vH1jrR1VMPdrh5XYuz61F4P/3HeTCaO4w7u0OGxRNLRNdGEZNO5fqLlZBlD8rIZkpfD0PxshuXnMDQ/h7FjhzJ8UC4jBuUycnAuIwfnMHJILqOG5DF6aB6jhuQxZlieHkcg8jbo/54kmXnKSL527YwjX7gQ+6J1wpdv3BcwxH0px9XBHSf25ezhy7ujw48cM/bF/eZxOvzN147wje/hy9w723c/Kgl0iv9qj08eR33l9/D973isx9cde+tbO5JU49azYokyy4ysLCPbYu+zs7LIzoptz8my2HszcrJjCTc7C3KyssjJNrKzjNysLHJzjJysLHKzs8jPCa+5sfW8nCwG52YzKDebQTnZDM7LJi9Ht6IRiYqSUJJMGTuUKe+aEnUYIiIpTf8EFBGRyCgJiYhIZJSEREQkMkpCIiISGSUhERGJjJKQiIhERklIREQioyQkIiKRsa63VpGjmVk1sCPqOBI0FtgXdRAnQKacB+hcUlGmnAek9rmc6u6FvVVSEsogZlbi7sVRx/F2Zcp5gM4lFWXKeUBmnIuG40REJDJKQiIiEhklocxyf9QBnCCZch6gc0lFmXIekAHnomtCIiISGfWEREQkMkpCacbMFpjZJjMrM7M7u9meb2a/DttXmNnk/o8yMQmcy8fNrNrM1oTlE1HE2Rsze9DMqsxsQw/bzczuDee5zszm9neMiUrgXC4xs7q4z+Rr/R1jIsxsopktN7ONZlZqZp/tpk5afC4JnktafC7dij2hU0s6LEA2sBWYCuQBa4EZXer8A/CTsH4j8Ouo434b5/Jx4IdRx5rAuVwMzAU29LD9GmApsQfLzgdWRB3z2ziXS4Ano44zgfMYD8wN68OBzd38/UqLzyXBc0mLz6W7RT2h9DIPKHP3be7eAjwGLOxSZyHwUFh/HLjczHp49nakEjmXtODuLwD7j1FlIfCwx7wMFJjZ+P6J7vgkcC5pwd0r3X1VWG8ANgJFXaqlxeeS4LmkLSWh9FIE7Ip7X85b/zIeqePubUAdMKZfojs+iZwLwHVhqORxM5vYP6GdcImea7q4wMzWmtlSM5sZdTC9CUPSc4AVXTal3edyjHOBNPtcOikJpZfuejRdpzcmUicVJBLnH4DJ7n4u8Gfe7OGlm3T5TBKxitjtWGYBPwB+F3E8x2Rmw4AngM+5e33Xzd3skrKfSy/nklafSzwlofRSDsT3BiYAFT3VMbMcYCSpObzS67m4e427N4e3PwXe0U+xnWiJfG5pwd3r3b0xrD8F5JrZ2IjD6paZ5RL70v6lu/+mmypp87n0di7p9Ll0pSSUXl4FppnZFDPLIzbxYEmXOkuARWH9euBZD1cuU0yv59JlfP59xMbC09ES4OYwG2s+UOfulVEH1RdmdnLnNUYzm0fsO6Qm2qjeKsT4ALDR3e/uoVpafC6JnEu6fC7dyYk6AEmcu7eZ2e3A08Rmlz3o7qVm9k2gxN2XEPvL+gszKyPWA7oxuoh7luC5/KOZvQ9oI3YuH48s4GMws0eJzU4aa2blwF1ALoC7/wR4ithMrDKgCbglmkh7l8C5XA982szagEPAjSn6j5wLgY8B681sTSj7CjAJ0u5zSeRc0uVzeQvdMUFERCKj4TgREYmMkpCIiERGSUhERCKjJCQiIpFREhIRkcgoCYkkiZk1vs39Hzezqb3Uec7Mit9unS71C83s/yVaX+TtUBISSUHh3l/Z7r6tv9t292qg0swu7O+2ZeBREhJJsvCL/P80sw1mtt7MPhzKs8zsx+EZMU+a2VNmdn3Y7aPA7+OOcZ+ZlYS63+ihnUYz+56ZrTKzZ8ysMG7zDWb2ipltNrOLQv3JZvaXUH+Vmb0zrv7vQgwiSaUkJJJ8HwRmA7OAK4D/DLck+iAwGTgH+ARwQdw+FwIr495/1d2LgXOBvzGzc7tpZyiwyt3nAs8Tu9tBpxx3nwd8Lq68Cnh3qP9h4N64+iXARcd/qiLHR7ftEUm+dwGPuns7sNfMngfOC+X/1907gD1mtjxun/FAddz7D5nZbcT+nx0PzADWdWmnA/h1WH8EiL/RZef6SmKJD2K34/mhmc0G2oEz4upXAacc53mKHDclIZHk6+mhgsd62OAhYBCAmU0Bvgic5+61Zvbzzm29iL8nV+fdyNt58//7zwN7ifXQsoDDcfUHhRhEkkrDcSLJ9wLwYTPLDtdpLgZeAV4k9tC+LDM7idiNQzttBE4P6yOAg0BdqHd1D+1kEbuRJcBHwvGPZSRQGXpiHyN2I9lOZwAbEjg3kbdFPSGR5Pstses9a4n1Tr7k7nvM7AngcmJf9puJPS2zLuzzR2JJ6c/uvtbMVgOlwDbgf3po5yAw08xWhuN8uJe4fgw8YWY3AMvD/p0uDTGIJJXuoi0SITMb5u6NZjaGWO/owpCgBhNLDBeGa0mJHKvR3YedoLheABa6e+2JOJ5IT9QTEonWk2ZWAOQB33L3PQDufsjM7gKKgJ39GVAYMrxbCUj6g3pCIiISGU1MEBGRyCgJiYhIZJSEREQkMkpCIiISGSUhERGJjJKQiIhE5v8DNjmXYBt85Z8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE on Training set : 753.7824867899045\n",
      "RMSE on Test set : 786.0172079616123\n",
      "r2_score on Training set : 0.8430423465498388\n",
      "r2_score on Test set : 0.8547372903751442\n"
     ]
    }
   ],
   "source": [
    "#1. 生成学习器实例，5折交叉验证\n",
    "lasso = LassoCV(cv = 5)\n",
    "\n",
    "#2.模型训练\n",
    "lasso.fit(X_train, y_train)\n",
    "alpha = lasso.alpha_\n",
    "print(\"Best alpha :\" , alpha)\n",
    "\n",
    "#3. 模型性能：cv\n",
    "mse_cv = np.mean(lasso.mse_path_, axis = 1)\n",
    "rmse_cv = np.sqrt(mse_cv)\n",
    "print(\"cv of rmse :\", min(rmse_cv))\n",
    "\n",
    "#4. 显示不同alpha对应的模型性能\n",
    "plt.plot(np.log10(lasso.alphas_), mse_cv) \n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show() \n",
    "\n",
    "#训练误差\n",
    "y_train_pred = lasso.predict(X_train)\n",
    "rmse_train = np.sqrt(mean_squared_error(y_train,y_train_pred))\n",
    "print(\"RMSE on Training set :\" ,rmse_train)\n",
    "\n",
    "#测试误差\n",
    "y_test_pred = lasso.predict(X_test)\n",
    "rmse_test = np.sqrt(mean_squared_error(y_test,y_test_pred))\n",
    "print(\"RMSE on Test set :\" ,rmse_test)\n",
    "\n",
    "r2_score_train = r2_score(y_train,y_train_pred)\n",
    "r2_score_test = r2_score(y_test,y_test_pred)\n",
    "print(\"r2_score on Training set :\" , r2_score_train)\n",
    "print(\"r2_score on Test set :\" , r2_score_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "season_1        -973.953630\n",
       "season_2          -0.000000\n",
       "season_3           0.000000\n",
       "season_4         474.247837\n",
       "mnth_1          -172.547973\n",
       "mnth_2           -68.043612\n",
       "mnth_3           264.912762\n",
       "mnth_4            -0.000000\n",
       "mnth_5           369.746549\n",
       "mnth_6             0.000000\n",
       "mnth_7          -400.514428\n",
       "mnth_8            15.141078\n",
       "mnth_9           643.203086\n",
       "mnth_10          332.445993\n",
       "mnth_11         -237.828881\n",
       "mnth_12         -179.219138\n",
       "weathersit_1     384.226728\n",
       "weathersit_2      -0.000000\n",
       "weathersit_3   -1431.529484\n",
       "weekday_0       -276.459777\n",
       "weekday_1       -134.851710\n",
       "weekday_2        -25.494016\n",
       "weekday_3         -0.000000\n",
       "weekday_4         49.312954\n",
       "weekday_5         78.658175\n",
       "weekday_6         84.137077\n",
       "temp            2942.057118\n",
       "atemp            891.718331\n",
       "hum            -1652.758660\n",
       "windspeed      -1398.540434\n",
       "holiday         -311.125900\n",
       "workingday         1.473680\n",
       "yr              1941.255795\n",
       "dtype: float64"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lasso模型增加了L1正则，系数值进行了收缩，同时有些特征的系数为0。由于岭回归模型在训练集上的误差和Lasso模型差不多，但在测试集上的误差比Lasso模型小。故岭回归模型比Lasso模型性能稍好。"
   ]
  }
 ],
 "metadata": {
  "_change_revision": 0,
  "_is_fork": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
